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ABSTRACT 

We study the post-main sequence stars in NGC 604, the most luminous HII region in 
M33. Previously, a number of Wolf-Rayet (WR) stars and one red supergiant (RSG) 
have been discovered. Based on broadband photometry of the region, we present ev- 
idence that is consistent with the presence of this RSG and with three more RSG 
candidates. Using SED fitting based on HST UVIJHK photometry we estimate the 
ages of the WR stars and RSGs finding that the two populations are from distinct for- 
mation episodes with ages 3.2 ±1.0 Myrs and 12.4 ± 2.1 Myrs, respectively. The RSGs 
have greater extinctions towards their line of sight than the WR stars consistent with 
the RSGs producing large amount of dust. Using the WR and RSG populations and 
similar SED fits to the most massive O stars we estimate that the total stellar mass 
is (3.8 ± 0.6) x 10 5 Mq. We find a large discrepancy between the expected Ha flux 
from such a massive cluster and that one observed. This suggests that 49^jg percent 
of the ionizing photons produced by massive stars in NGC 604 is leaking from this HII 
region. We also suggest that the implications of an old RSG population mean that if 
NGC 604 was more distant and only observed in the infrared (IR) it would be difficult 
to study the youngest burst of star formation due to the contamination of RSGs. 

Key words: stars: evolution - binaries: general - stars: Wolf-Rayet - stars: super- 
giants - HII regions: NGC 604 - galaxies: individual: M33 



1 INTRODUCTION 

NGC 604 is the most luminous and massive HII region in 
the nearby galaxy M33. Due to its large size, great mass 
and proximity it has been intensively studied to gain in- 
sight into the physics of HII regions and the evolution 



of m as sive stars (e.g . Hunter et alj 1 19961: Terlevich et al 



1996: Drissen et al.l 



E. mum 
19931. 



2008 



Mafz-Apellaniz et al 



2004 : iGonzalez-Delgado fc Perez! |200d . among others'). Of 



paramount importance is to determine the age of the stellar 
population and the different stellar types present in this re- 
gion, because a deeper knowledge of this object can help us 
to infer properties for star-forming regions of similar types 
in distant galaxies. 

The confirmed stellar constituents of NGC 604 are 
OB stars, WR stars and at least one spectroscopi- 
cally confirmed RSG dHunter et al.1 1 19961 : iTerlevich et all 
1 19961 : iDrissen et ail 120081) . The different stellar types that 
NGC 604 harbours makes this HII region an ideal object 
to study its stellar population and to test stellar evolu- 
tion and population synthesis models. NGC 604 is a key 
region to test whether the population inferred from re- 
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solved stellar populations can also match the unresolved 
spectral features that would be used in more distant 
galaxies. In this paper we use this region to qualitatively 
check the accuracy of the binary pop ulation and spectral 
synth esis (BPASS) code described in lEldridge fc Stanwavl 
(2009) and at http://www.bpass.org.uk. Here we ana- 
lyze the resolved stellar population observed with deep 
Hubble Space Telescope (HST) images and study the re- 
gion a s an unresolved object by using the integrated spec- 
trum ( Gonzalcz-Dcl gado fc Perezll2000! ) and observations of 
iRelafio fc Kennicuttl (|2009h . 

The possible presence of RS Gs in NGC 604 come s 
from the spe c trosco pic evidence of ITerlevich et all l|l996l ). 
iHunter et al.l (|l996l ) suggested the existence of a supergiant 
population based on HST stellar photometry, but these au- 
thors could not elaborate further these ideas since their ob- 
servations wh*reJumted_tothe visible part of the spectrum. 
Recently. llBarba et al.1 (2009) have shown the complexity of 
the stellar population of this region using infrared NICMOS 
(Near Infrared Camera and Mult i- Object Spectrometer) ob- 
servations. They found evidence of massive young stellar ob- 
ject candidates and suggested the existence of a RSG pop- 
ulation within the region. 

Here we aim to study in detail the post-main sequence 
objects RSGs and WR stars of NGC 604, along with the 
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most massive main- sequence stars in the region. We ob- 
tain the photometry for the stars within NGC 604 in the 
IR-NICMOS filters and com bine the results with the opti- 
cal photometry performed by [Hunter et all 1 19961 ). We ana- 
lyze the Spectral Energy Distribution (SED) of these stars 
from the optical to the IR wavelength range and investi- 
gate whether the RSG and WR populations have different 
ages. We then attempt to infer from these observations the 
total stellar mass and check if this is consistent with other 
observations such as the total Ha flux from NGC 604. 



2 EVIDENCE OF RED SUPERGIANTS IN 
NGC 604 

The stellar content of NGC 604 was previously studied by 
Hunter et al. (1996) using UVI photometry. These authors 
showed the existence of numerous massive main-sequence 
stars and identified a candidate WR population in their 
colour-magnitude diagram (CMD). They also suggest the 
existence of numerous stars bright enough to be blue or red 
supergiants, which would indicate an older sub-population in 
NGC 604. Due to the limitation of their data, restricted only 
to the visible part of the spectrum, they did not study this 
cool and evolved stellar population. Taking advantage of the 
available N1CMOS-IR data for this region and the compari- 
son with the optical photometry, we now investigate further 
this evolved population. Using Hunter et al.'s photometry 
we selected stars in the centra l cluste r of the region (Cluster 
A, as defined in lHunter et al. (1996)) with magnitudes and 
colours in the range of F555W-F814WM and F555W<-6 
and make a list of possible RSG candidates in NGC 604. 
These limits ensure we are taking massive (M > 15 Mq) 
RSGs and therefore we can compare this popula tion with 
the massive WRs detected bv lDrissen et all (|2008l 1. 

2.1 NICMOS data analysis 

NGC 604 was observed with the NICMOS (NIC2) on board 
the Hubble Space Telescope (HST) under the ID program 
10419 (IP: R. Barba) with the filters F110W, F160W and 
F205W (similar to J, H and K filters). We retrieved 6 fields 
from the HST Multimission Archive at the Space Telescope 
Institute (MAST), 5 covering the central part of NGC 604 
and another covering an outside field for background sub- 
traction. Each image has a field of view of 19.2 x 19.2 arcsec 
and a pixel scale of 0.075 arcsec. 

The retrieved data were processed by the pipeline soft- 
ware Calnica (Version 4.4.0) but further analysis needed to 
be performed to obtain the final mosaic images. First we, 
masked out the Coronagraphic Hole and applied the Mul- 
tiDrizzle software packageQto create the final cosmic ray- free 
mosaic images. The retrieved F205W images showed some 
residuals due to a bad flat-field correction after being pro- 
cessed by Calnica with the default calibration files provided 
in the data archive. Therefore we generated a new flat-field 
using the F205W background images and reprocessed the 
images with it. The use of this flat-field clearly improved 
the final images for this filter. Finally, the ~20 rows of the 



CCD affected by vignetting were eliminated from the im- 
ages. We then created the mosaic image for the 5 fields in 
each filter using the World Coordinate System information 
from the headers. The final image for the F205W filter is 
shown in Fig.[T] The angular resolution (FWHM) in the fi- 
nal images are 0.16, 0.18, 0.20 arcsec for the corresponding 
filters F110W, F160W, F205W. 

We obtained the stellar magnitudes for each filter us- 
ing t he crowded-field fitting software DAOPHOT (|Stetsonl 
1 19871 ) . We identified stars with intensities above 5 times the 
r.m.s. noise of each image using DAOFIND, then aperture 
photometry was performed on the selected stars with a circu- 
lar aperture of 3 pixel radius (0.23 arcsec). A Point-spread 
function (PSF) was fitted using bright and isolated stars 
in each image. Then, the instrumental magnitude was de- 
rived for each identified star using the fitted PSF in each 
image. We subtracted the identified stars off the original 
image and we run DAOFIND again in the residual image 
to select fainter stars that were not identified in the origi- 
nal images. The magnitudes of these new stars were derived 
and included in the final photometry. In this second itera- 
tion of DAOPHOT we detected the following stars in each 
filter, 22 stars in F205W, 41 stars in F160W and 47 stars 
in F110W. However, we are only modelling stars detected 
in the 3 NICMOS filters. With this condition the second it- 
eration add only 5 stars to the sample in the models. None 
of these 5 stars are in the group classified as WR, RSG or 
massive OB stars, and therefore it does not have impact in 
the analysis nor conclusions presented in the paper. Aper- 
ture corrections to a radius of 0.49 arcsec were applied to 
the instrumental magnitudes in order to use the photomet- 
ric calibrations recorded in the image headers and provided 
by the NICMOS teanfl Finally, magnitudes in the Vega 
system for each star and each filter were then obtained. 

The result of the stellar photometry is shown in the 
top panels of Fig.[2] We present here the colour-colour and 
CMD for the identified stars in NGC 604. The magnitude 
errors are a combination of the image photon noise and the 
uncertainties derived from the PSF fitting procedure and 
from aperture corrections; these are shown as a function of 
the magnitude for each filter in the lower panels of Fig. [2] 
The red diamonds in the top panels of Fig. [2] correspond to 
the RSG candidates derived from the UVI photometry from 
iHunter et al.l |j~996;). We identified the RSG candidates in 
the NICMOS images using the coordinates of Hunter et al.' 
star catalogue and the F555W image from their study. Only 
6 candidates were within the field of view of NICMOS. With 
exception of one RSG candidate (number 2 in Table[2]), for 
which we identified another star very close to it with mag- 
nitudes in all NICMOS filters and with no counter part in 
the UVI photometry, the rest of the stars are clearly identi- 
fied as single stars in the NICMOS images. For the star 2 in 
Table[2]we added the IR fluxes of the two identified stars in 
order to perform the SED fitting with all the filters available 
(see section 3). 
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Figure 1. Final mosaic F205W image of NGC 604 obtained with NICMOS-NIC2 Camera. The total field of view of NICMOS observations 
is 50"x60" with a pixel scale of 0.075 arcsec. 



2.2 Confusion with foreground stars 

The identification of RSGs in extragalactic clusters is a diffi- 
cult task as they could be confused with foregroun d Galactic 
stars. As it has been explained in iMassevl |l998j, a 40 M 
RSG in M33 would have an apparent magnitude in V-Band 
between +15.5 to +18, depending on its effective temper- 
ature. For dwarf late-type stars the absolute V magnitude 
is +6 (for a KO dwarf star) and +12 (for a M5), taking 
~ 600 pc as the estimated distance for these stars in the 
direction of M33, the apparent V magnitudes for these stars 
will be between +15 and +21 i|Massevlll9"9Sh . Therefore, it 
is possible to mistake a Galactic dwarf late-type star for a 
RSG in M33. 

A method of separating foreground stars from RSGs 
based on IR two-colour diagrams have been proposed by 



lElias fc Frogel (|l985f ): the RSG are 0.1-0.2 mag redder than 
late-type dwarfs in J-H at comparable H-K. In Fig. [3] we 
show the JHK colour-colour diagram for the stars with NIC- 
MOS photometry as well as for the RSG candidates. The 
magnitudes in the NICMOS niters have been converted into 
magnit udes in the JHK syst em following the transforma- 
tion of iBrandner et aL I i|200ll ). Since the magnitude conver- 
sion procedure already adds uncertainties in the magnitude 
we only convert those stars with magnitudes errors lower 
than 0.2 mag. The blue line denotes the location of the fore- 
ground late- type dwarf stars with magnitudes obtained from 
iKoornneeH (| 1983h . Clearly four of our RSG candidates stand 
between 0.2 to 0.6 magnitudes above the location denned by 
the foreground stars. The two stars that lie below the line 
defined by the foreground stars are S2 and S5. These stars, 
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Figure 2. Top panel: Colour-colour diagram (left) and CMDs (centre and right) for NGC 604 in the NI CMOS filters F110W , F160W 
and F205W. The red diamonds are the candidate RSG stars selected using the UVI stellar photometry from lHunter et ahl l ll996l ). Bottom 
panel: Observed magnitude as a function of the estimated errors for each filter: F110W (left), F160W (centre) and F205W (right). No 
extinction corrections have been applied to the magnitudes and all the magnitudes are in the Vega system. 



as we will see later, have anomalous SED fitting (see sec- 
tion 4). 



3 STELLAR EVOLUTION MODELS AND SED 
FITTING 

Determining the physical parameters of stars from pho- 
tometry requires a combination of stellar evolution models, 
model stellar atmospheres and a method to find the best 
fit between models and theory. In this section we outline 
the stellar and atmosphere models, their amalgamation and 
how we use the results to estimate parameters for the ob- 
served stars. To calculate the absolute magnitudes of the 
obse rved stars we take the distance of NGC 604 to be 840 
kpc (|Freedman et al.lll99ll ). 

We use stellar models from the Cambridge STARS code 
(Eldridge, Izza rd fc Toutll2008l . and references therein) that 
have been correlated with stellar atmosphere models to pro- 
duce synthetic spectra and broad-band pho tometric colours 
as described in lEldridge fc Stanwavl l|2009l ) . Their key fea- 
ture is that there is not only a set of detailed single star mod- 
els, but also an extensive set of detailed binary star models 



which are key to producing a realistic synthetic stellar pop- 
ulation. We consider stell ar models at the relevant metallic- 
ity fo r NGC 604 Z=0.008 jVflchez et al.lll98g| ; lEsteban et all 
|2009|) . where Z is the metallicity mass fraction (a metallic- 
ity of Z=0.020 is conventionally considered solar). It is well 
known that local metallicities within disk galaxies can vary 
by significant factors away from the mean host metallicity. 
We inspected the affect of varying metallicity by running the 
fits with models at metallicities of Z = 0.004 and 0.020. We 
find that our results change only slightly and allowing the 
metallicity to vary would only increase the uncertainty of our 
results. We assume a hydrogen mass fraction of X = 0.73 
and a helium mass fraction of Y = 0.262. We calculate 
observable magnitudes for these stellar interior models by 
using a model atmosphere that best represents the appear- 
ance of the interior model at each timestep. The atmosphere 
models are for OB stars, WR stars and other stellar type s 
and they are described in Smith. Norris fc Crowtherl J2002D 
Ham ann. Grafener fc Liermannl |200rj ) and We stera et al.l 
|2002l ). respectively. 

Given that stellar evolution is non-linear, and binary 
evolution is even more so, we do not interpolate between 
models with different initial parameters to determine the 



The RSGs & WR stars of NGC 604 5 




Figure 3. Colour-colour diagram in the JHK system for the stars 
identified in NICMOS images and magnitude errors less than 
0.2 mag. Red diamonds correspond to the RSG candidates shown 
in Fig. [2] The blue line shows the location of foreground late- 
type dwarf stars. The transformation from the NICMOS filter to 
J HK system has been p erformed following the equations given 
in iBrandner et al.l 11200111. The reddening vector corr esponds to 
Av=5mag and the ICardelli. Clayton fc Mathisl dl989h 
law. 



extinction 



nature of the observed stars. Instead, we compare the ob- 
served SEDs to the synthetic SEDs at every timestep of each 
stellar model and estimate a most likely mass, age and ex- 
tinction for each observed star. This is done by calculating 
mean parameters from all the models that match the ob- 
served SED. 

We first estimate the amount of extinction required for 
the synthetic SED to match the observed SED. We calculate 
the amount of extinction required in each band to achieve 
an exact match between the observed and model magni- 
tudes. Using the ICardelli. Clayton fc Mathisl (| 19891 ) extinc- 
tion curve we convert each magnitude to the equivalent ex- 
tinction in A v calculating the conversion factor , such that 
Ai = ctiA v for the ith filter. We then calculate the mean ex- 
tinction from the six filters, as follows, 



1 

A v ,„(m,t) = 



(M 



i{m,t)) 



(1) 



where the magnitudes predicted by a model are 
M m odei,i(m,t) for a star with initial mass m and age t, the 
observed magnitude is M b s ,n,i for the nth star in the ith 
filter and N is the number of filters in which the source is 
detected. This mean extinction is then used to correct the 
synthetic SED magnitudes for dust. 

This dust corrected model magnitudes, 

M m (m,t) = M modc i,i{m,t) +aiA VtTl (m,t), (2) 

are then used to calculate the probability, p n (m,t), of a 
match between the observed SED and theoretical SED can 
be evaluated as follows: 

6 



Pn(m,t) = JJexp 



(■^modcl,cor,i,n (^i ^) -^^obs,n,i) 
^ u total,i 



(3) 



calculated as a combination in quadrature of the error in the 
photometry and an assumed error in the stellar models of 
0.3 mag. This error is estimated from the magnitude differ- 
ence between two stars of neighbouring masses in our grid of 
models. We must also take into consideration the resolution 
of our mass grid in the probability calculation. For example 
our 15-Mq model has neighbouring model masses of 14 and 
16M©. While the 80Mq model has neighbours in the mass 
grid of 70 and IOOMq. Therefore we add a factor, Am(m), 
to the probability to give every mass an equal probability. 
For the the 15 and 80Mq examples this factor would be 1 
and 15 respectively. 

The last factor to consider in the likelihood of a match 
is how quickly stars evolve. Two stars might have the same 
apparent SED but one is evolving more rapidly than an- 
other. A star evolving slowly, such as main-sequence stars 
will be more probable matches to a SED that a short-lived 
WR stars. Therefore we use the probability and the timestep 
of each model SED as a weight to estimate values for the pa- 
rameters of each observed star. For example, the mass (M n ), 
age (r n ) and extinction (-A v , n ) are obtained as follows: 



M n = 



Em Et " m At ( m > t ) Am ( m ) Pnjm, t) 
El' Et " At(?n, t) Am(m) p n (m, t) 

Ejv Et" t At(m, t) Am(m) p n (m, t) 
Em Et 11 At (m, t) Am(m) p„ (m, t) 



(4) 



(5) 



(7) 



The errors are therefore assumed to be Gaussian. Ototal i s 



A _ El' Ef ^y,n ( m > 1 ) At(m, t) Am(m) p„ (m, t) 
Em Et 11 At (m, t) Am(m) p n (m, t) 

and for each of these we are able to estimate an error from 
the variance, for example: 

^ 2 = Em Et" m2 At ( m ' f ) Am(w)p„(m, t) _ m2 
Em Et" At(m, t) Am(m) p„(m, t) 

Thus for each observed SED we obtain an estimate of age, 
initial mass, extinction, luminosity, surface temperature, the 
mass ratio (q — M2/M1) for the binary models and an error 
for each value. We note that we are only determining three 
independent parameters from each fit: luminosity, surface 
temperature and extinction. The age and mass are model 
dependent and linked to the surface temperature and lumi- 
nosity. If we re-ran the analysis with different or updated 
models the luminosity, surface temperature and extinction 
would remain similar but the age and mass would change 
by a greater degree. 

We perform our fitting process using two sets of mod- 
els: single star models and binary models, then we com- 
pare the probability from the two fits to determine which 
is more likely. We find, except for the RSG1, 4 and 6, the 
binary models give a better match for all the stars we fit. 
This is likely to be due to the fact the binary models can 
cover more of the Hertzsprung-Russel diagram than the sin- 
gle stars therefore have a greater chance of achieving a close 
match to the SED. In most cases the likelihood difference 
between single and binary models is minimal up to only fac- 
tors between 1 to 3. The only two fits for which binaries are 
significantly favoured are RSG3 and WR1. For these two 
stars they are 700 and 30 times more likely to be binaries 
respectively. 

The SED fits for the RSGs and WR stars assume that 
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the stars are post-main sequence objects and therefore we 
only use post-main sequence models to prevent the fits being 
affected by main-sequence models skewing the fit. This is 
because in weighing by the timestep a WR star and a main- 
sequence OB star might have very similar SEDs but a star 
is more likely to be a main-sequence star due to its slow 
evolutionary timescale. 

In Figure [4] we show examples of the location of possi- 
ble SED matches on the Hertzsprung-Russel diagram. Then 
Figure [5] compares the theoretical SEDs from the same fits 
to the observed SED. For the red supergiant the surface 
temperature is tightly constrained but the luminosity is less 
well constrained. The more luminous models require more 
extinction to achieve the same SED and so tend towards 
higher surface temperatures. However, despite this large ap- 
parent degeneracy, a low luminosity around 10 L is pre- 
ferred because of the consideration of the timestep. While 
there are many possible values of the luminosity, there is 
only one place where the star may be observed for an ap- 
preciable period of time. Therefore the fitting process gives 
greatest weight to this region and therefore gives the masses 
and ages listed in Table [2] below. In this Figure we also in- 
clude a similar plot with the timestep weight removed from 
the SED fit. Without the timestep weight a better match 
is achieved by the more luminous and hence more massive 
RSG. 

We note that in Figure [4] we predict more mas- 
sive RSGs than other stellar evolution models such as 
iMevnet fc Maederl ||2003h . However the lifetimes of these 
massive RSGs are very short due to their high mass-loss 
rate and so they are unlikely to be observed. Furthermore 
the evolution of these massive RSGs is also uncertain and 
extra p ulsation driven mass-loss may shorten their lifetimes 
further l|Yoon fc Cantiello|[2010l ). 

For the WR star the region on the surface temperature 
is less well constrained. This is not surprising as these stars 
are hot and with the filters used in the SED fitting we are 
only looking mainly at the Rayleigh- Jeans tail of the black- 
body spectrum so there is little constraint on the surface 
temperature. Therefore the SED fits of the WR and main- 
sequence stars are intrinsically less well constrained. Fur- 
thermore in Figure [S] the agreement between the observed 
and theoretical SEDs is slightly poorer, especially in the J 
and H filters. This is due to most of the WR flux being 
output at shorter wavelengths so the IR flux may be more 
greatly effected by IR emission from other nearby sources. 

An important check to understand the accuracy of the 
SED fitting process is to feed theoretical models through the 
SED fitting code and compare the output to known inputs. 
The models we use to do this where calculated separately 
to have initial masses between the mass grid of our stel- 
lar models. We use single star models with initial masses of 
65, 90 and 1OOM . We also include a 15OM model to un- 
derstand how stars above our maximum model mass would 
appear. We calculate fits to the models at three evolutionary 
phases, the main-sequence, RSG phase and the WR phase. 
The results of these fits are shown in Table [T] 

The fits on the RSG phase are most accurate as it is 
possible to achieve a better match for cooler stars. This is 
because the peak of black-body spectrum is observed rather 
than for hotter stars where only the Rayleigh- Jean tail of 
the black-body spectrum is covered. Thus both the main- 



Table 1. Results of test of SED fits for stellar models not in our 
standard grid. The evolutionary phases are, MS main sequence, 
RSG - red supergiant and WR - Wolf-Rayet. 



Model 


Single 


Binary 




mass 


star fit 


fit 


Evolutionary 


/M 


M,/M 


Mi/M @ 


Phase 


65 


85 ±26 


79 ±25 


MS 


65 


60 ±5 


55 ± 9 


RSG 


65 


97 ±21 


83 ±29 


WR 


90 


63 ± 15 


54 ±17 


MS 


90 


98 ± 18 


103 ± 19 


RSG 


90 


90 ±20 


85 ±28 


WR 


110 


89 ± 25 


98 ±24 


MS 


110 


88 ± 14 


89 ±21 


RSG 


110 


85 ± 25 


80 ±31 


WR 


150 


87 ±24 


85 ±25 


MS 


150 


106 ± 10 


104 ± 17 


RSG 


150 


85 ± 25 


80 ±31 


WR 



sequence and WR fits are less accurate. However it is pos- 
sible to explain the inaccuracy of the the 65 and 8OM fits. 
The the 65M SED was taken from towards the end of evo- 
lution for the star while the 9OM SED was taken from 
the near beginning of the track. Therefore there is some de- 
generacy in the fitting of main-sequence stars between mass 
and age. In addition the more massive stars, including the 
15OM do tend to have their masses underestimated due to 
this mass-age degeneracy. Only with more information from 
the UV of the stars could we achieve a more accurate fit. 
An important conclusion here is that while individual fits 
will be uncertain by using the properties of all of the stars 
some of these errors will cancel out to provide a reasonable 
estimate of the mass and age of NGC 604. 

This test indicates that the information inferred on the 
RSGs can be trusted more than that on the hotter stars. 
Also that while each individual fit may not give the exact 
result, if we take the entire population of stars into consid- 
eration our combined uncertainty in the total age and mass 
of a region is significantly reduced. 

Another problem with our fit is due to weighing by the 
timestep a star in a rapid phase of evolution may be misiden- 
tified as a lower mass, more slowly evolving, star. Therefore 
we removed the At from the fitting process. We found that 
the resulting values were similar to those when the timestep 
is included. The only exception was RSG1. In Table[2] we 
show that the mass of the RSG changes from 23 ± 20 M to 
67 ± 36 M with a corresponding drop in age and increase 
in extinction. Such a massive RSG has a very short evolu- 
tionary timescale which is why it is not favoured when the 
timestep weight is included. It is not possible from SED fit- 
ting alone to determine which mass is the best solution even 
though the more massive, shorter lifetime, RSG fit is less 
probable. To determine its true identity a spectrum of the 
star is required. However we assume its has a similar age to 
the other RSGs in our analysis of the region and we note 
that its location (Fig. [9} is close to the centre of NGC 604 
and suggests that the more massive solution is possible. 
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Figure 4. Examples of the region on the theoretical Hertzsprung-Russell diagram where models match the observed SEDs of RSG1 
and WR7. We include two plots for RSG1 showing the two fits with and without timestep weighing, RSGla and RSGlb respectively. 
The dashed grey lines indicate single star stellar models with masses of 10, 20, 30, 40, 50, 60, 70, 80, 100 and 120Mq. The grey points 
indicate the regions where a SED match is made, the darker the point the better the agreement between the observed and model SED. 
The square with green error bars indicates the SED fit with the highest value of mAt(t)p„ (m, t) from both the single star and binary 
models. The triangle with red error bars indicates the mean T e ff and L estimated from the single models and the diamond with blue 
error bars indicates the mean T cf f and L estimated from binary models. 



RSGla 



RSGlb 



WR7 






Figure 5. Comparison of the observed and theoretical SEDs for the two solutions of RSG1 and WR7. The dashed grey lines are 
theoretical matches to the observed SED, the darker the line the higher the probability of a match. These lines are the SEDs for the grey 
points in Figure[4] The green line is the theoretical SED with the highest value of mAt(t)p n (m, t). The red line indicates the mean SED 
estimated from the single models and the blue line indicates the mean SED estimated from binary models. These three lines correspond 
to the colours points in Figure(4] The asterisks represent the observed SED, the vertical lines represent the uncertainty in the observed 
magnitudes. 



4 RED SUPERGIANTS IN NGC 604 

From Fig.[2]we find that there are at least four RSGs within 
NGC 604. Determining their ages is vital to understand- 
ing the star-formation history of this region. The results of 
the SED fitting are shown in Table[2] We see that, ignor- 
ing the unlikely younger fit for RSG1, the RSGs all have 
ages in the range of 10 to 15 Myrs with masses between 
14 to 23 Mq. Th ese are the typical a ges and masses ex- 
pected for RSGs (|Levesciue et al.l l2005). One important re- 
sult is that RSGs 1, 4 and 6 have a large extinction excess 
above that one predicted fr om the extinction map derived in 
iRelano fc Kennicutd l|2009h . In Fig.0 we plot the expected 
extinctions of the stars versus the values derived from SED 
fitting. Each of these stars is luminous with a bolometric 



magnitude in the F205W filter (similar to K-band) of ~ —10. 
Such extinction exce ss is found to be ty pical for RSGs with 
similar luminosities (|Massev et al.ll2005h . The excess is ex- 
plained as intrinsic extinction from dust produced in the 
RSG atmosphere and stellar wind. A final detail to note is 
that RSG3 is significantly better fit by a binary model than 
a single star. 

The second solution for RSG1 with a younger age and 
higher mass also has a greater extinction excess. Again this is 
not unusu al when compared to the excesses found in Galac- 
tic RSGs (|Massev et al.l [2005V However, if this solution is 
true, then this star could be the most luminous RSG in 
the local group, as recent observations indicate that the 
maximum observ ed luminosity for ty pical RSGs is around 
log(L/L ) « 5.3 jMassev et all 120051 ). Although due to its 
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Table 2. Results of the SED fitting of the most massive RSG and WR stars in NGC 604. We also show the results for the most luminous 
stars in the Hunter's catalogue. Note that RSG2 and RSG5 may also be foreground stars and we list the two possible solutions to the 
SED fit for RSG1. 
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higher surface temperature it would probably be classified as 
a yellow supergiant (YSG). Its mass and age would be con- 
sistent with the hotter stellar population described below. 
However, only a spectrum of the object will truly confirm its 
nature, as it is the case for RSG4, the star sp ectroscopically 
confirmed as a RSG in lTerlevich et all ljl99rJ l. 

The two impostors RSGs S2 and S5 are most likely to 
be foreground stars as discussed above. The SED fit gives 
masses for the stars similar to the WR and OB stars in 
NGC 604. Their large extinctions are similar to the magni- 
tude of the RSGs so they may be transition objects either 
becoming WR stars or leaving the main-sequence to become 
the RSGs. Only a spectrum for these objects can reveal their 
true nature. 

In summary, the mean age of the RSGs 1, 3, 4 and 6 
is 12.4 ± 2.1 Myrs and the mean initial mass is 18 ± 4Mq. 
The other two suspected RSGs are not RSGs but probably 
foreground stars. 
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5 WOLF-RAYET STARS IN NGC 604 



In NGC 604 there are a number of s pectroscopically con- 
firmed WR stars (jDrissen et al.l [2008) . From the available 
photometry we are able to perform a SED fit for 10 of these 
objects and the results are shown in Table[2] The derived 
ages and masses are remarkably consistent, with only WRs 
5 and 8 being significantly different from the mean mass and 



Figure 6. The age and mass values derived from the SED fitting 
for the RSGs (red), WR stars (green) and main-sequence stars 
(blue). The lighter blue points are based on UVI SED fitting 
alone. The names of the RSGs and WR stars are given. 
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Figure 7. Extinction comparison between the values obtained 
from the Balmer extinction map of iRelano fc KennicutH ||2009| ) 
and those derived her from the SED fitting. The observations are 
colour codes so that the RSGs are red, WR stars are green and 
main-sequence stars are blue. The names of the RSGs and WR 
stars are given. 




Figure 8. The theoretical Hertzsprung-Russcll diagram for 
NGC 604. The single-star models in grey have masses of 20, 30, 40, 
60, 80, 100, 120 M . Colour code is: red for RSGs, green for WR 
stars and blue for main-sequence stars. The lighter blue points 
are based on UVI SED fitting alone, names for the RSGs and 
WR stars are given. The Rib star is the second more luminous 
fit for RSG1. It also has a slightly higher surface temperature. Its 
location on the Hertzsprung gap means if this is a solution the 
stars is in a rapid phase of evolution. 



age of all the WR stars. For each of the stars there is also 
a best fit with massive companions with a mass ratio of 0.5 
to 0.8. 

The WR stars have little excess extinction. We plot the 
expected extinction of the stars versus the values obtained 
from SED fitting in Fig. [7] We see th at the WR stars cor- 
respon d with the extinction derived in lRelano fc Kennicuttl 
(2009). This is to be expected as any dust around the star 
created during its brief RSG phase will be quickly swept 
away by the fast and dense WR winds (e.g. Eldridge ct al. 
2006). The case with the highest extinction, WR6, is iden- 
tified as a WNL star. There stars are thought to have still 
some of their remaining hydrogen envelope and therefore 
will not have completely removed all the dust from their 
immediate vicinity. This is consistent with the evolutionary 
scenario of a star being a RSG before becoming a WNL star. 

The consistent ages and masses of the fits to the WR 
stars suggest that they were all formed in a single and recent 
burst of star-formation with a mean age of 4.0 ± l.OMyrs. 
This is significantly different to the RSG ages shown above 
but it is also similar to the age inferred for RSG2 and RSG5 
and the alternative fit for RSGT The mean initial mass for 
the WR stars is 76 ± 10M Q . 



6 MAIN-SEQUENCE STARS IN NGC 604 

In the photometric catalogue of iHunter et al.l (|l996l ) there 
are 1040 stars with UVI photometry. We have performed a 
similar process of SED fitting for the entire catalogue using 
these three filters. Because most of these stars are main- 
sequence objects the inclusion of JHK photometry does not 
increase the accuracy of our fit. However, for a sample of the 
most luminous stars in the catalogue we have also performed 
a 6 filter fit and listed the results in Tabled The estimated 
parameters agree within the fit errors. From the fit of the en- 
tire catalogue, we find 34 O stars with masses above 70 Mq , 
with the most massive star » 100 Mq . For these 34 stars we 
find a typical age of 2.4 ± 0.3 Myrs. Below 70 Mq the error 
in the age of the SED fit increases significantly, as shown in 
Fig.O because of the degeneracy between mass and age for 
O stars. For example above 60 Mq we find 53 stars with a 
mean age of 2.8 ± TlMyrs. 

In most cases of the OB stars the derived extinction is 
si milar to that seen across th e NGC 604 region as derived 
bv lRelano fc Kennicuttl ([2009). This is as expected with the 
OB star winds being fast and sweeping out the region around 
each star. 

In both mass range mentioned above the ages are 
slightly younger than the ages derived for the WR stars but 
still overlapping within the errors. This is to be expected be- 
cause both the OB and WR stars come from the same initial 
mass range but are respectively earlier/later in their relative 
lifetimes. We created synthetic cluster models of NGC 604 
formed in a bursts lasting between 1 to 2 Myrs. We find that 
in these models the true age of the cluster is between the 
ages of the OB and WR star populations with the same ini- 
tial masses. In addition the duration of the burst is similar 
to the spread of the two ages of the two populations. 

It is worth noting that the apparent spread also has 
another possible explanation the stars might have different 
rotation rates with those rotating more rapidly having longer 
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main-sequence lifetimes. If this is the case then the duration 
of the star-formation burst would be decreased. 



7 STELLAR POPULATION OF NGC 604 

From the analysis shown above we see that NGC 604 is a 
collection of main-sequence stars, RSGs, possible BSGs and 
WR stars. We find two distinct populations: an older one 
with an age of 12.4 ± 2.1Myrs represented by the RSGs, 
and a younger one which is a combination of O stars, pos- 
sible massive RSG/BSGs and WR stars with a mean age 
of 2.4 ± 0.3Myrs for the O stars and 4.0 ± l.OMyrs for the 
WR stars. We take the age of this younger population to be 
3.2 ± l.OMyrs by combining these last two measurements 
assuming that the true age of the burst is intermediate be- 
tween the relatively young O and old WR stars. Importantly 
the WR stars indicate the upper age for the star formation 
burst to be around 5Myrs. This can be visually checked in 
Figure HJ] with the ages clustered closely around 2 to 4 Myrs. 

We plot our results for the estimated luminosities and 
surface temperatures onto a Hertzsprung-Russell (HR) dia- 
gram in Fig. [5] We see that many of the WR stars cluster 
around the track for the 80 Mq , along with many of the most 
massive stars and S2 and S5. While RSG 1,3, 4 and 6 are 
closer to the 20 Mq track. This again demonstrates that the 
RSGs in NGC 604 are a separate population from the more 
massive stars. 

For the young and old populations it is possible to place 
limits on the total stellar mass formed in each population. 
Using the mean mass and its uncertainty to describe the 
mass range of the stars in each population, we can estimate 
the total mass of the stars assuming a Salpeter IMF (with 
a — —2.35) from 0.5 to 1 20 Mq and a slope with a — —1.3 
between 0.1 and 0.5Mq (|Kroupall2002l ). This only gives us 
the mass of primary stars, therefore we multiply the resul- 
tant mass by 1.5 to account for each star having a binary 
companion. Doing this we assume a flat distribution of the 
mass ratio between the two stars with < M1/M2 ^ 1. 
For the older population with 4 RSGs assuming they rep- 
resent all stars between 14 and 22 Mq we find a total mass 
of 1700 ± 900 Mq. For the younger population with masses 
above 70 M , we estimate a mass of (3.8 ± O.6)xlO 5 M . 
We note that if we use all stars with masses above 60 Mq 
we obtain an identical mass estimate. 

The stellar population of NGC 604 can also be investi- 
gated by using various spectral f eatures in a integrated spec- 
trum for all stars in the region. iGonzalez-Deleado fc Perea 
l|2000l) obtained UV and optical spectra of the NGC 604 
region and studied in detail the nebular emission of the re- 
gion. Using their spectrum we have calculated the equiva- 
lent widths for the stellar population diagnostic lines of CIV 
at 1500A, Hell at 1640A and the Blue WR bump around 
4686A and compared them to those predicted by our syn- 
thetic spectra for the stellar populations. The calculation 
of the equivalent widths and the con struction of the syn- 
thetic spectra are described in detail in lEldridee fc Stanwavl 
(|2009T) . 

We plot the comparisons in Fig. 1101 We see that the 
features that are dependent on the WR stars -the Hell line 
and Blue WR bump- are in agreement with the observed 
values are ages around 3 Myrs. The CIV line which is dom- 



inated by emission from the main-sequence OB stars has 
agreement at much earlier ages of 1.5 Myrs. Therefore we 
can assume that there is some agreement between models 
and the observed spectra but there is no exact match. We 
suggest this is because the two spectra do not fully cover 
the whole area of NGC 604 and therefore we cannot be sure 
that the entire resolved population of stars is contributing 
to the spectrum observed. Also the model spectra are for an 
instantaneous burst where NGC 604 has an extended period 
of star-formation. However, the general agreement provides 
confidence in the synthetic binary population and spectral 
models. 



8 Ha AND LEAKAGE OF IONIZING 
PHOTONS 

Using the masses for the stellar population given above and 
the predicted Ha per Mq at different ages listed in Table[3]it 
is possible to estimate the Ha flux these two populations will 
produce. We calcula te the nebula em ission from Cloudy 
Ferland et alj [l998) as described in Eld ridge fc Stanwavl 



20091 ). The ionizing spectrum given to Cloudy as input 



is obtained using the fitted stellar populations in this pa- 
per (Table[2]) and the corresponding stellar interior and at- 
mosphere models. The predicted flux varies most with the 
assumed age for the populations. For the old population 
we find log(_F(Ha)/ergs _1 ) ^ 36.75, while for the younger 
population of 3.2 ± l.OMyrs we find log(F(Ha)/ergs~ 1 ) = 
39.921q2J. Therefore, it is the younger population the one 
which dominates the Ha emission from NGC 604. 

The Ha flux from NG C 604 has been studied in detail by 
I Relano fc Kennicuttl (|2009h . They obtained an extinction- 
corrected Ha of log(F(Ha)/ergs _1 ) = 39.63 ± 0.07. Our 
model prediction agrees within the errors with this value, 
but in general overestimates the Ha flux and we investigate 
this further. 

iRelafio fc~K cnnicutt (2009), using the single-star Star- 
burst99 code ( Leitherer et al.l Il999l ) inferred the mass of 
stars in NGC 604 to be 5 x 10 s Mq assuming a Salpeter 
IMF between 0.1 and 120Mq. This value was derived using 
single star models and an age of 4 Myrs. Using the observed 
flux, our single star models and an age of 4 Myrs we also 
find a similar mass of 4.7 x 10 5 Mq. However, the uncer- 
tainty in the mass is greater for single star models because 
for them the Ha flux varies more with age than for binary 
star models (see Tabled . The weaker variation of Ha flux 
for binaries compared to that from single star models is due 
to rejuvenation of secondary stars and mergers from binary 
evolution. 

The binary population models produce more Ha flux 
per stellar mass over a longer period of time. As most of 
the stars in Table[2]are better fit with the binary population 
models we favour the Ha flux prediction from these models. 
The difference in the predicted and observed Ha flux gives 
us an estimate of the leakage fraction of ionizing photons. 
For NGC 604 we derive a leakage fraction of 491JJ per cent 
of the total ionizing flux in the region. 

This value of leakage is dependent on uncertainties 
within the stellar models. Alternative explanations for the 
mismatch of observed and predicted Ha flux are therefore 
that either we are overestimating the number of massive 
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Figure 9. Final mosaic F110W image of NGC 604 with the WR stellar population (blue circles) from lHunter et all jl996l ^ and the RSG 
stars (red circles) identified in NICMOS data. 



stars or the age of the region must be greater than 5Myrs. 
Both these possibilities reduces the Ha flux from the stel- 
lar population. However, the large number of young main- 
sequence stars and massive WR stars indicates such an old 
age is unlikely. Therefore we can be confident in our results. 



9 DISCUSSION & CONCLUSIONS 

We have found further evidence that there is a popula- 
tion of RSGs within NGC 604. This older population in- 
ferred from the RSGs 12.4 ± 2.1Myrs. The age of the 
more numerous young stars is 3.2 ± l.OMyrs, which is 
in agreement with previous studies (iMaiz-Apellaniz et al.l 
l2004l ; lGonzalez-Delgado fc PeredlioOCh . 

There is some uncertainty in the nature of RSG1 as it 



may be a very massive RSG or YSG in a short-lived phase 
of evolution. To accurately determine its mass a spectrum of 
this object is required. However its deep position in NGC 604 
surrounded by many O and WR stars indicates that it may 
well be a very massive RSG. 

The RSGs have larger extinctions than the other stars 
in NGC 604. This is in agreement with them producing 
large amounts of dust to obtain greater intrinsic extinc- 
tion (iMassev et al.l 120051 ). The one WNL WR star also 
has a higher extinct i on th an expected from the study of 
iRelano fc Kennicuttl (|2009l ). This is in agreement with the 
evolutionary scenario that RSGs lose their hydrogen en- 
velopes and become WR stars. The two stars labeled in Ta- 
ble[2] as S2? and S5? may be RSGs in the very late phases 
of losing their hydrogen envelopes due to their very high 
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Figure 10. The equivalent widths of different stellar population diagnostic lines of NGC 604, the dashed lines with errors represented 
by the hashed region, compared to single-star (black) and binary (red) population models. Left: the CIV line at 1550 A from O-stars, 
centre: the Hell line at 1640 A from WR stars and right: the Blue WR bump around 4686 A. 



Table 3. Predicted Ha flux per Mq of stars for both single star 
and binary populations. 
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inferred extinction. However they are more likely to be fore- 
ground stars. 

We also find that the WR spectral features Hell and 
the Blue Bump and the O-star CIV P-Cygni profile of the 
composite spectrum of the region agree in general with the 
results from the resolved stellar population. The OB main- 
sequence stars appear to be younger than the WR stars. 
However, a direct comparison is difficult as the UV and op- 
tical spectra of the region do not cover the entire region and 
the duration of the star-formation burst is similar to the age 
of the cluster and therefore difficult to model. 

By using the estimated age of the stellar population we 
are able to predict the Ha flux from the massive stars and 
compare this to the flux observed for NGC 604. We find 
the two agree in magnitude but that 49lj| percent of the 
ionizing flux may be escaping from the region. The great- 
est uncertainty in calculating the Ha flux, and therefore 
the amount of leakage, is the age of the stellar population. 
Thus, we are unable to completely rule out that there are 
no absolute leakage of ionizing photons. However, since leak- 
age has been stu died to be a common property of luminous 
HII regions (e.g. lOev fc KennicuttJ[l997l : IZurita et al.ll2002l; 
iRelano fc Peimbertj20o3 T it is more likely that leakage is oc- 



curring. Future study of the main-sequence population using 
HST/ACS data will enable tighter constraints to be placed 
upon the predicted ionizing flux escaping from NGC 604 
(Laporte et al., in prep). 

The evidence supporting the existence of RSGs in 
NGC 604 and their relation to the most recently formed 
stars is important because recent and future cutting-edge 
telescopes operate at these longer wavelengths (e.g. Spitzer, 
Herschel, JWST and the E-ELT). Soon all premier obser- 
vatories will operate predominantly at IR wavelengths and 
therefore modelling and observation efforts must act now 
to address this observational shift. T here are known prob- 
lems with modelling near-IR emission ( Maraston et al . 2006; 
iLanonl [20081 ): typically it is found that for some star clus- 
ters ages derived from near-IR observations do not agree 
with those estimated from optical ones. The RSGs within 
NGC 604 are a key element required to study distance HII 
regions in the IR. The relative contributions of older RSGs 
and younger stars must be evaluated. Key to breaking any 
degeneracy will be the use of spectra and spectral lines such 
as the IR lines of RSGs a nd Wolf-Rayet stars (e.g. lCrowtherl 
l2007l : lDavies et al.ll20TOl ). This implies that the study of dis- 
tant HII regions at IR wavelengths will be challenging unless 
the contribution of such contaminants can be fully under- 
stood. 
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